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ABSTRACT 



Aims. Neutron stars (NSs) produced in the Milky Way are supposedly ten to the eighth - ten to the ninth, of which only ~ 2 x 10' are 
observed. Constraining the phase space distribution of NSs may help to characterize the yet undetected population of stellar remnants. 
Methods. We perform Monte Carlo simulations of NS orbits, under different assumptions concerning the Galactic potential and the 
distribution of progenitors and birth velocities. We study the resulting phase space distributions, focusing on the statistical properties 
of the NS populations in the disk and in the solar neighbourhood. 

Results. It is shown that ~ 80 percent of NSs are in bound orbits. The fraction of NSs located in a disk of radius 20 kpc and width 
0.4 kpc is < 20 percent. Therefore the majority of NSs populate the halo. Fits for the surface density of the disk, the distribution of 
heights on the Galactic plane and the velocity distribution of the disk, are given. We also provide sky maps of the projected number 
density in heliocentric Galactic coordinates (/, b). Our results are compared with previous ones reported in the literature. 
Conclusions. Obvious applications of our modelling are in the revisiting of accretion luminosities of old isolated NSs, the issue of 
the observability of the nearest NS and the NS optical depth for microlensing events. These will be the scope of further studies. 
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Neutron stars are born during core-collapse of massive (M > 
8Mo) stars (type lb, Ic and II supernovae, herafter SNe) or, less 
frequently, by accretion-induced collapse of white dwarves that 
reached the Chandrasekhar's limit (type la SNe). A rough es- 
timate of the total number of NSs generated in the Milky Way 
(MW) can be obtained from the present-day core- collapse SN 
rate, B ^s , which is of the order of a few per century dPiehl et aLl 
120061) and assuming (i^s constant during the lifetime of the 
Galaxy (~ 10 Gyr). Hence the total number of NSs created in 
the MW lies between 10^ and 10''. NSs may then represent a non 
negligible fraction of the Galactic stellar content. 

Up to now only ~ 2 x 10^ NSs have been observed, the ma- 
jority of which as isolated radio pulsars (PSRs) with ages far 
shorter (< 100 Myr) than the MW hfetime. Older NSs have 
been detected only when recycled in binary systems by mass 
and angular momentum transfer from a companion star, thus be- 
coming millisecond pulsars (see Lorimer 2008 and references 
therein). Isolated old NSs (ONSs) have not been identified so 
far because, once their energy reservoir, both thermal and rota- 
tional is exhausted, they are pretty close to being invisible. As 
a consequence little is known about their physical and statistical 
properties. 

On the other hand the expected phase-space distribution 
of ONSs can be constrained by means of population synthe- 
sis models once a realistic set of initial conditions is given. 
Population synthesis studies of G alactic NS have been per- 
formed by many authors in the past. iHartmann et al.l (1 19901) and 



iPackzvnskil (1 19901 hereafter H90 and P90 respectively) stud- 
ied the orbits of Galactic NSs, looking for a possible link with 
gamma-ray bursts. These studies differed in their assumptions. 
In particular H90 assumed a Gaussian distribution centered at 
200 kms ' for the distributi on of NS birth ve locities, while P90 
adopted the distribution of L vne et al] (Il982h which gives a dif- 
ferent (higher) weight to the low velocity tail of the distribution. 

NS orbits in the Galactic gravitational potential were also 
investigated by Blaes andRajagopal (1991), Blaes and Mada^ 
(T993), Zane et al,, (119951) . iPopov and Prokhorovl 01998 ) and 
Popov et al. (2000, hereafter BR, BM, Z95, PP98 and POO re- 
spectively) in order to constrain the number of nearby NSs ac- 
creting from the interstellar medium (ISM, see' Treves et alj20^ 
and references therein). BR, BM and Z95 adopted initial condi- 
tions similar to those of P9 except th e distribution of birth ve- 
locities which, following Naravan and Ostriker (1990) , was as- 
sumed to be Maxwellian with a dispersion of 60 km s~' . 

'Popov et all (l2000h explored the observability of accret- 
ing ONSs for a wide range of initial mean velocities, be- 
tween and 550 kms"', assuming a Maxwellian distribution. 
The paucity of observed accr etorjj in the ROSAT catalogue 
jNeuhauser and Trumpet] IT999I) led to the conclusion that NSs 
are born with average velocities of at least 200 km s"' . 

This is confirmed by observations of known young NSs. 
PSRs show in fact spatial velocities of several hundreds kms ', 
i.e. of the same order of the escape velocit y from t he MW (see 
e.g. lLorimer et aljfl 99711 Arzoumanian et al . 2002., Brisken et al.1 



No sources have been positively identified so far. 
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I2003L iHobbs et alj [20051 and iFaucher-Giguere and Kaspill2006t 
hereafter L97, A02, B03, H05 and F06 respectively). Some PSRs 
exhibit velocities in excess of 1000 km s"' . A striking example is 
PSR B 1508 +55: the proper motion and parallax measurements 
obtained from radi o observations points to a transverse velocity 
of - 1083 km s ' (IChatteriee et"aLll2005h . 

Similar high values of the velocity have been inferred also 
for objects belonging to other classes of isola ted NSs. Thanks to 
Chandra observations. lHui and Beckerl(l2006l) estimated a veloc- 
ity of ~ 1100 kms~' for the centra l compact object RX J0822- 
4300. Recentlv lMotch et al.l (l2009l) measured the proper motion 
of one of the ROSAT radio-quiet, thermally emitting NSs (the 
Magnificent Seven) and found a value of the 3D velocity of 
600 - lOOOkms"'. This is not uncommon in PSRs and hence 
they concluded that the velocity distribution of the Magnificent 
Seven is not statistically different from that of normal radio pul- 
sars. 

The origin of such high velocities is not at all clear. 
An asymmetric SN explosion is c onsidered one poss i ble ex - 
planation (e.g. 'Shklovskiij 1 19701 iDewev and Cordes' 1 19871) 
Also the effects of binary disruption (e.g. Blaaw Il961l 
llben and Tutukovl IT996I) may contribute to the observed ve- 
locities. Recently it has been proposed that the fastest NSs 
are the remnants of runaway progenitors expelled via N- 
body interactions from the dense core of yo ung star clusters 
dGvaramadze. Gualandris and P orte gies Zwarll 2008). 

If aU classes of isolated NSs share the same typical birth ve- 
locities, no matter how these are achieved, a large fraction of 
these objects can escape the potential well of the MW in a rela- 
tively short time. This fact has consequeces of all observable NS 
populations. In this paper we focus on the effect of high birth 
velocity and likely evaporation on the still elusive population of 
ONS. Constraining the expected phase space distribution is in 
fact crucial to define suitable strategies for thier detection. 

Based on recent estimates of the birth velocity distribution, 
in this paper we reconsider the dynamics of isolated NSs. We 
perform integration of stellar orbits using our new code PSYCO 
(Population SYnthesis of Compact Objects), developed for this 
purpose. In Section |2] we describe the ingredients of the simu- 
lation, i.e. the gravitational potential of the Milky Way and the 
distributions of progenitors and birth velocities. We present the 
results of the simulation in Section [3] In particular we investi- 
gate the statistical properties of the NS population in the Galactic 
disk and at the solar circle. We fit the surface density of the disk 
and the average height distribution and compute the surface and 
volume densities in the solar vicinity. We fit also the velocity 
distribution in the disk, both with respect to the Galactic center 
and the local rest frame of the ISM. We discuss our results and 
their possible implications in Section|4] The results of this work 
will constitute the base for further studies on the observability of 
ONSs. 



2. METHOD 

We follow an approach similar to P90. Initial conditions (posi- 
tion, velocity) are taken randomly from the selected distributions 
and assigned to each synthetic NS by means of a Monte Carlo 
procedure. 



2.1. Distribution of progenitors 

The initial positions of NSs in the Galaxy are defined in a galac- 
tocentric cylindrical coordinates system {R,<p,z), where the z 



axis corresponds to the axis of rotation of the MW. These ini- 
tial positions reflect the distribution of NS progenitors: accord- 
ing to Bronfan et al. ( 2000, hereafter BOO), formation of massive 
stars is currently concentrated in a annular region which follows 
the distribution of molecular hydrogen. However, to explore the 
effects of different initial conditions on the current phase-space 
configuration of NSs, we choose four possible radial distribu- 
tions of progenitors from the literature. 

P90 adopted an exponential probability distribution, based 
of the observed sur face brightnesfl of face-on Sc galaxies 
(Ivander Kruitlll987h 



R { R \ 
p(R) dR = fls-Tj— exp I - - — dR , 



(1) 



where p{R)dR is the probability that a NS is born between R and 



R + dR, Re 



4.5kpcandaR = 1.0683. 



BOO obtained the already mentioned radial distribution of 
starforming regions in the MW from the combined far infrared 
(FIR) and millimetric emission produced by clusters of massive 
stars embedded in ultra-compact HII regions. The FIR (surface) 
luminosity, p{R), has a Gaussian shaped rise until it reaches a 
maximum at ~ 4.7 kpc (FWHM of 2.38 kpc), then it decays 
exponentially, with a scale-length of 1.78 kpc, towards the outer 
part of the MW, approximately following the distribution of neu- 
tral hydrogen. The radial birth probability is then obtained from 
the FIR luminosity from the equation 



p(R) dR = 



Rp{R)dR 
Rp(R)dR 



(2) 



Another possible distribution of NS progenitors can be 
obtained from the surface d ensity of Galactic SN remnants 
(ICase and BattacharyaLll998L hereafter CB98) 



exp 



m - Ro) 



(3) 



where a - 2 and y6 = 3.53 and Rq - 8.5 kpc. The corresponding 
radial probability density is again obtained from equation (|2]l. 
The fourth radial distribution adopted has been proposed by 

F06 



V27ro- 



exp 



(R - Rpeak) 



(4) 



where Rpeak - 7.04 kpc and cr - 1.83 kpc. This distribution 
ha s been extrapolated from th e observed PSR distribution found 
bv Yusif ov and Kucukl(l2004l) . Finally for all models we assume 
that NSs can be born from to 20 kpc. 

It is our opinion that the distribution proposed by P90, in 
spite of being obtained from observations of external galaxies, 
may better represent the long term star formation history of the 
MW. The other models are based on the present-day distribution 
of population I objects, which could have been rather different in 
past epochs (see for example Chiappini et at. 2001). The models 
of BOO, CB90 and F06 are probably better suited for population 
studies of young/middle-aged NSs (PSRs, magnetars etc.). 



In th J band. 
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Fig. 1. Normalized radial probability distribution of NSs pro- 
genitors. P90 (solid line), BOO (dotted), CB98 (dashed) and F06 
(dot-dashed). 

Table 1. Parameters of the spiral arms. 



Arm 


k 


R. 


00 






[kpc] 


[radians] 


Norma 


4.25 


3.48 


1.57 


Carina-Sagittarius 


4.25 


3.48 


4.71 


Perseus 


4.89 


4.90 


4.09 


Crux-Scutum 


4.89 


4.90 


0.95 



2.1.1. Spiral arms and initial distribution of heights 

Massive stars are located in the spiral arms of the MW (they are 
indeed the ideal tracers of the spiral structure), thus we model 
spiral arms in the distribution of NS progenitors adopting the 
same prescription of F06, i.e., NS progenitors are distributed 
along four logarithmic spirals, each spiral described by the equa- 
tion 

<I>{R) = k MR/R,) + </>o . (5) 

The values of the parameters k, Rt and (po for each spiral are 
given in Table [1] Actually, equation |5] describes the position of 
arm centroids. A more realistic distribution can be obtained if 
the positions of progenitors are scattered, both in the radial and 
azimuthal directions, around these centroids (Fig.|2|. Details on 
how the scatter is added to the initial positions of NSs can be 
found in F06. 

Th e thickness of the star forming region is few t ens of parsecs 
(B00.[Maiz-Apellanizll2001h . However, as P90 and lSun and Hani 
(l200 4l) have pointed out, the long term dynamical behavior of a 
NS populatio n is insensitive to the scale height of its progeni- 
tors (see also Kiel and Hurlevll2009h . Following their results we 
assume that all NSs are born on the Galactic plane (z - 0). 

2.2. Distribution of birth veiocities 

The true form of the distribution of birth velocities is still a hotly 
debated issue and few constrains exist (see F06 for an exhaustive 
discussion). HP97, L97 and H05 proposed a Maxwellian distri- 
bution 





-70 1- I- I I i 

-20 -10 10 20 

X [kpc] 

Fig. 2. Initial positions of NSs - Radial distribution from P90. 
The position of the Sun is (8.5, 0.0). 



Alternatively, iFrver et al.1 (Il998l) . ICordes and Chernofj (Il998l) . 
A02 and B03 proposed a bimodal distribution 




where w is the relative weight of the two components of the dis- 
tribution. 

Using the same PSR sample of B03, F06 explored, together 
with Maxwellian and bimodal models, other possible distribu- 
tion functions like the double-sided exponential 

Mv,)= y^exp(-^), (8) 

where v,- represents a single component of the velocity and Vexp 
is a characteristic velocity; the Lorentzian 

/'(V,) = i r , (9) 

ny[\ + (vf/r^)) 

where y is a scale parameter defining the half- width at half m ax- 
imum, and the distribution proposed by iLvne et all ( Il982h and 
adopted by P90 

Mv)=— (10) 

where again v represents the tridimensional velocity. F06 con- 
cluded that the Maxwellian model is less favored. On the other 
hand they disfavour also the bimodal distribution and prefere in- 
stead single parameter models, pointing out that the bimodality 
found by other authors arises if thse alternative single parameter 
models are not investigated. 

To explore the effects of the birth velocities on the fi- 
nal phase-space distribution of NSs, we adopt the Maxwellian 
model of H05 as well as four of the models proposed by F06, 
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Fig. 3. Differential velocity distributions obtained from sim- 
ulated velocity vectors. H05 (solid), F06DG (dotted), F06E 
(dashed), F06L (dot-dashed) and F06P (triple dot-dashed). 

Table 2. Velocity distribution models. 



Model 


Parameters 


< v> 






[kms"'] 


H05 


cr = 265 km s"' 


420 


F06B 


tri = 160 km s"' 


335 




0-2 = 780 km s"' 






w = 0.9 




F06E 


v„p = ISOkms"' 


400 


F06L 


7 = 100 km s"' 


447 


F06P 


V. = 560 km s"' 


331 



i.e. the bimodal, the double-sided exponential, the lorentzian and 
that of P90. From here on we refer to these models as H05, F06B, 
F06E, F06L and F06P respectively. 

The value of the mean tridimensional velocity for each dis- 
tribution is calculated numerically from simulated velocity vec- 
tors (Table |2]i. All the velocity distributions refer to the Local 
Standard of Rest (LSR) of NS progenitors. Thus, the true 3D ve- 
locity of a neutron star with respect to the Galactic Reference 
Frame (GRF) is the vector sum of the birth velocity and the cir- 
cular velocity at the birthplace, v - y birth + Vc,>c- 

2.3. Gravitational potential 

Once the initial conditions have been assigned, the motion of 
NSs is described by the equation 



-V<D, 



(11) 



where r = r(/?, 0, z) is the position of the NS and <l> is the grav- 
itational potential of the MW. We adopt the same 3-component 
model of Smith et al. ( 2007. hereafter S07) 



(12) 



where Og, O/j and O^ represent the bulge, disk and halo contri- 
butions respectively. 

The gravitational potential of the bulge is jHernquistl 1 1 990h 



B 



(13) 



where Mb = 1.6 x 10'" and rg - 0.6 kpc are respectively 
the mass and the scale radius of the bulge and r - y/R^ + 7} is 
the distance from the GC. 



The disk potential ha s instead the following form 
dMivamoto and Nagailll975h 



On = -- 



(14) 



,2 



where Mo =5x10 Mq is the mass of the disk and the - \ 
kpc and zu - 0.3 kpc are respectively the scale length and scale 
height of the disk. 

Finally, the potential of the halo is 
(iNavarro. Frenk and Whit3.ll996l) 



47r Gp,r„. /, cr\ 



(15) 



where 



Ps = 



ln(l -Hc)-c/(l H-c) 



(16) 



is the characteristic density, c is the concentration parameter, r,,,> 
is the virial radius and pcr is the critical density of the Universe. 

The parameters of potential are the same of S07, except for 
the concentration parameter c and the virial radius r,,,> (19.2 and 
274 kpc respectively), which were adjusted to match the lAU 
standard values for the distance of the Sun from the Galactic 
center, /?() = 8.5 kpc, and the circular velocity at the solar circle, 
i'c/>c(^o) = 220 km s"', together with the escape velocity from 
the MW at the same distance, VesciP-o) - 544 km s"' (S07). The 
corresponding value of the virial mass. 

Very recently ^ Reid et al] (12009 ) gave a new estimate of the 
circular velocity, v„>e(^o) - 254km s"' with - 8.4 kpc, 
meaning that the MW may be more massive that previously 
thought. To asses the effect of the enhanced mass of the Galaxy 
on NS orbits, we choose a further set of parameters for the po- 
tential: the masses of the bulge and disk are increased by a fac- 
tor (254/220)^, i.e. the ratio of the squared circular velocities in 
the two cases. For the halo, the concentration parameter c re- 
mains the same while the virial radius r,.„. is in this case 332 
kpc, which yields an 80 percent increase of the virial mass, 
M,„> ~ 1.8 X IO'^Mq. 

We note that in the model where Vc„e(^()) = 254 km s ' we 
get Vesc - 664 km s"' . This is higher than the central value (544 
km/s) estimated by S07; however, it is not far from their 90% 
upper limit (608 kms"'), especially when we consider that Vg^c 
was obtained by assuming v'cve = 220km s"', and that modify- 
ing such assumption introduces further uncertainty in its deter- 
mination (Smith, private communication). 



2.4. Orbit Integration 

We calculate the orbits of 150000 NSs for each model (Table 
O, assuming that NSs are born at constant rate during the whole 
MW lifetime (10 Gyr). Hence the age of each NS is selected 
randomly from an uniform distribution. The orbit of each NS 
is then calculated via numerical integration of the system of 
equations ( fTTT ). for a time corresponding to its assigned age. The 
axial symmetry of the potential implies conservation of angular 
momentum with respect to the axis of rotation of the MW. This 
allows to reduce the number of equations in (fTTI ) to four 
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Fig. 4. Rotation curve for our Milky Way model (solid). Dotted, 
dashed and dot-dashed represent the bulge, disk and halo contri- 
butions respectively. 
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Fig. 5. Escape velocity on the Galactic plane. The circular veloc- 
ity (dashed) is plotted for comparison 



dR 
dt 
dz 
dt 
dvR 
dt 
dvz 
dt 



-r = vr. 



50 



(17) 



where j, is the angular momentum with respect to the z 
axis. Integration of equations ( fTTll is performed wi th a 4th or- 
der Runge-Kutta algorithm (e.g. iPress et aTl 1 19921) with cus- 
tomized adaptive stepsize. The relative accuracy of integrations 
is kept below 10"^ using the energy integral E as reference, i.e. 
i6E/E) < IQ-^, where 



£ = _ + (D(r). 



(18) 



To limit the computation tim^ and avoid lockups of the 
code, all NSs reaching 0.1 kpc from the Galactic center are dis- 
carded. The fraction of NSs traveling to within 0.1 kpc from the 
Galactic center is less than 1 percent in each run. 



Table 3. Models for initial conditions. (*) denotes models with 
updated potential. 







Birth velocity distr. 




Spatial 


H05 


F06B 


F06E 


F06L 


F06P 


distr. 












P90 


lA 


IB 


IC 


ID 


IE 


BOO 


2A 


2B 


2C 


2D 


2E 


CB98 


3A 


3B 


3C 


3D 


3E 


F06 


4A 


4B 


4C 


4D 


4E 


P90 


lA* 


IB* 


IC* 


ID* 


IE* 
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0.001 
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15 
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Fig. 6. Surface density of NSs in the disk, obtained from best fit 
parameters (A^^^,. = 10^). Models IB (soHd line), 2B (dotted), 
3B (dashed) and 4B (dot-dashed). 

3. RESULTS 

Our calculations show that the statistical properties of NSs are 
affected mostly by the distribution of birth velocities while the 
effects of different distributions of progenitors are less promi- 
nent. For this reason we focus on results of models lA to IE, i.e. 
with the distribution of progenitors of P90. Results of models 
differing only for the distribution of progenitors are quite sim- 
ilar, the only substantial difference is the shape of the surface 
density (Fig. in fact, in models based on the P90 progenitor 
distribution, the density peaks at the center (Rpeak - 0) whereas 
for other models the density peaks off center (2 < Rpeak ^ 6 
kpc). 

3.1. Fraction of bound neutron stars 

We first compute the fraction of NSs in bound orbits, /hound- 
Neglecting all those processes that could alter its energetic state 
(e.g. two body interactions), the final fate of a NS is known once 
its initial position and velocity are fixed. A NS star is bound 
when its initial velocity is lower than the escape velocity at the 
birthplace, v, < V(.5c(r), with 



v,.„.(r,) = V-2<l>(r,) , 

where r, is the position of the newborn NS. Thus 

N(V < Vesc) 



fboimd — 



(19) 



(20) 



The CPU time for a typical run is about 1 day. 



The retention fraction is ~ 0.7 for models lA, IC and ID, 
while for models IB and IE, fbound ~ 0.9 and 0.8 respectively 
(Table©. 
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Fig. 7. Distribution of heights /fzj - Model lA. Dashed Hne rep- 
resents the fitting function. 



3.2. Distribution of heights 

From here on our resuhs are obtained rescaling Nstar from 
150000 to lO**, which is our reference value for the total number 
of NSs produced in the MW. 

We study the distribution of NSs, /(z), as a function of the 
height on the Galactic plane. We adopt a logistic function 



/(z) = 



1 



{bob] + b2) ' 



(21) 



as fitting function (see e.g. Fig. |7]l. From these fits we estimate 
the average half density half thickness zi/2 of the disk (Table 
m. The values of the coefficients of the fit for each model, to- 
gether with the corresponding maximum error, are listed in the 
Appendix (Table |A.2t . The half density half thickness shows 
substantial variations from model to model, going from 100 to 
~ 400 pc for models lA to ID: for model IE in particular, z\/2 
is ~ 30 pc, i.e. roughly an order of magnitude smaller than those 
obtained from other models. We will return on this fact later. 
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Fig. 8. Final distribution of NSs in the disk - Model IE*. 
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3.3. Neutron stars in the disk 

Here we define the Galactic disk as the cylindrical volume with 
radius 20 kpc and height 0.4 kpc (i.e R < 20 kpc and |z| < 0.2 
kpc). The fraction of NSs that reside in the disk, fdisk, goes from 
~ 0.05 to ~ 0.20. Hence the majority of NSs born in the MW 
populate the halo (Tabled. 

We fit the logarithmic surface density of the disk adopting a 
fourth order polynomial as fitting function 

log I.(R) = flo + ajR + a2R^ + UjR^ + a4R'* . (22) 

The accuracy of these fits is always better than 5 percent. The 
values of the coefficients aj are listed in the Appendix (Table 
\AA\ . 

We made a visual check of the final distribution of NSs in the 
disk, looking for traces of the spiral arms. We found no evidence 
of spiral structure in the evolved distribution (compare Figs. [8] 
and|2]i. 

The hypothesis of constant formation rate yields an average 
age of ~ 5 Gyr, with young/middle-aged (< 10 Myr) NSs rep- 
resenting only ~ 0.1 percent of the whole population. As ex- 
pected this value is higher inside the disk, by a factor ~ 10, since 
NSs are born there, and did not have enough time to run away. 
However the excess of young NSs (Fig. |9]) does not alter the 
mean age in the disk, which is also ~ 5 Gyr. 



Fig. 9. Distribution of ages - Model lA. Solid an dashed lines 
represent global and disk populations respectively. 



3.3.1. Mean velocities 

The mean velocity of NSs in the disk is roughly the same for 
all models, < v > ~ 210 - 230 km s ' in the GRF while in the 
LSR the mean velocity is lower, < v"* > ~ 150 - 190 km s"' . 
An exception are models based on the distribution F06P, which 
show mean velocities in the LSR of ~ 80kms This fact can 
be easily explained: in the F06P model low birth velocities have 
higher probability (see Fig. O and thus the main contributor to 
the velocity of the star is the circular velocity, \hirth+'^circ - Vc,>c- 
The low velocity in the LSR implies also that NSs can not 
move too far away from the disk and that is the reason why, for 
models F06P, the scale height is considerably lower than in other 
models. 

Following Z95, we fit the cumulative velocity distribution 
of NSs, both with respect to the GRF and the LSR, with the 
following function 

(vIvS 

Giv) ^ \ \„ , (23) 

1 -(- V/V'o 
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Table 4. Statistical properties of NSs in the disk. 



Model 


lA 


IB 


IC 


ID 


IE 


lA* 


IB* 


IC* 


ID* 


IE* 


f bound 


0.70 


0.88 


0.72 


0.71 


0.79 


0.84 


0.91 


0.82 


0.77 


0.85 


fdisk 


0.05 


0.11 


0.10 


0.12 


0.19 


0.06 


0.13 


0.12 


0.16 


0.21 


Z\p 

[pc] 


367 


225 


164 


100 


33 


345 


192 


149 


80 


28 


<v> 
[km s"'] 


230 


220 


215 


213 


213 


262 


250 


249 


245 


245 


[km s"'] 


180 


146 


199 


164 


82 


199 


156 


216 


176 


89 


/so 


0.003 


0.001 


0.003 


0.002 


« 0.001 


0.002 


< 0.001 


0.002 


< 0.001 


« 0.001 


ho 


0.019 


0.060 


0.033 


0.074 


0.437 


0.014 


0.039 


0.029 


0.068 


0.379 



where a vq is a characteristic velocity. Fit values for vq, m and n 
are listed in the Appendix (Tables IA3] and lA!4b . 

3.3.2. The solar neig[nborhood 

To compare our results with previous works we focus now on the 
statistical properties of NSs in the so-called solar region, 7.5 < 
R < 9.5 kpc, as in BR, BM and Z95. 

From the fits we discussed in the previous subsection, we 
obtain the local surface density Zo which varies from ~ 0.4 - 2 x 
10^ (A^5,fl,./10')kpc"^. The volume density in the solar vicinity, 
no, also varies by a factor 5 between models, from ~ 1 - 5 x 10"'' 

We can now infer the distance of the nearest NS simply by 
computing the minimum volume around the Sun which contains 
at least a NS, assuming constant density 

1 = -r-«Ot/^ ^t/,m« ^\-, ; (24) 

3 \47rno/ 
typical values of c/„„„ are around 10 pc. 

3.4. Halo neutron stars 

The distribution of NSs in the halo is shown in Fig. [TO] Bound 
NSs can be found as far as ~ 1.5 Mpc, however the majority lies 
within the virial radius of the MW (~ 270 kpc). The radial dis- 
tribution clearly shows that unbound NSs start to be dominant at 
~ 500 kpc. Accordingly the mean velocity, after an initial drop, 
starts to rise almost linearly from ~ 500 kpc. The gravitational 
effects of other galaxies (e.g. Andromeda and the Magellanic 
Clouds) have not been considered. 

3.5. Density maps 

We compute the projected number density of NSs in heliocen- 
tric coordinates (I, b, d) and give the relative sky maps for stars 
within 30, 10 and 3 kpc respectively. Our sky maps (Fig. [T2] l 
clearly show that the most promising direction to look for NSs is 
towards the Galactic center, where the density is higher. Moving 
away from the center, the density drops rapidly even along the 
Galactic plane (FigfTTt. 



10° 




[ . . , , . . , , . . . ] 

2 4 6 

Radiol Distance [Mpc] 



Fig. 10. Model 1 A - (upper panel) Radial distribution of NSs in 
the halo, (lower panel) Mean velocity. 



In Table |6] we list the inferred values of the projected den- 
sity towards specific lines of sight (LOS). Tha sampling distance 
varies according to the LOS: for example, towards the GC the 
sampling distance, d„,ax, is equal to Ro while for the 3 other LOS 
lying on the Galactic plane c/„,„i is equal to the distance at which 
the LOS itself crosses the outer border of the stellar disk (20 
kpc). 

For large values of d„ax the projected density has non- 
negligible values even at high Galactic latitudes. One intriguing 
consequence is that NSs in the halo may contribute to the ob- 
servable rate of microlensing events, both of Galactic and extra- 
galactic sources (Sartore et al., in preparation). We thus calculate 
also the expected number density of NSs in the direction of the 
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Table 5. Statistical properties of NSs at the solar circle. 



Model lA IB IC ID IE lA* IB* IC* ID* IE* 

So 0.44 1.02 0.90 1.21 1.93 0.63 1.31 1.17 1.55 2.33 

(^)pc-^] 



no 1.1 2.6 2.3 3.0 4.8 1.6 3.3 2.9 3.9 

109 



rf™„ 12.9 9.8 10.2 9.2 7.8 11.5 9.0 9.3 8.5 7.4 
[pc] 

<v> 216 213 204 206 216 248 242 234 238 249 
kms"'] 

cv"''^> 173 140 191 158 72 191 150 203 167 79 
kms"'] 

/so 0.004 0.001 0.003 0.001 « 0.001 0.002 < 0.001 0.002 < 0.001 « 0.001 

/j'-ff'^ 0.025 0.043 0.035 0.084 0.497 0.014 0.032 0.034 0.078 0.454 



Line of sight Density [(^^) deg "] 

lA IB IC ID IE lA* IB* IC* ID* IE* 

/ = 0°, fc = 0° 0.42 0.90 0.84 1.01 1.08 0.61 1.16 1.05 1.22 1.26 
[X 10''] 

= 90°, Z^ = 0° 0.80 1.64 1.22 1.51 2.40 1.18 2.05 1.53 1.85 2.76 
[X 10'] 

= 180°, ^> = 0° 1.74 3.29 2.72 3.46 5.12 2.49 4.30 3.52 4.07 5.96 
[X 10^] 

= 270°, ^> = 0° 0.82 1.62 1.22 1.53 2.37 1.16 2.04 1.51 1.89 2.75 
[X 10^] 

Z7=+90° 5.12 2.83 2.86 3.55 3.63 0.76 2.84 4.39 6.47 1.48 
[X 10] 

fo = -90° 2.13 4.25 4.26 5.73 2.18 4.55 4.97 1.46 2.16 2.96 

[X 10] 

280°, b = -33° 3.25 3.35 3.47 2.87 2.46 3.79 3.39 4.32 3.31 2.66 

[X 10^] 

303°, b = -44° 4.44 4.25 3.76 2.78 2.78 5.20 3.85 4.07 3.02 2.79 



[10^ I 
[10-* 



Table 6. Projected density of NSs towars different lines of sight. 



/ 

I- 

/: 

/ = 

/ = 

[X 10^] 

IVIagellanic Clouds, assuming 48 and 61 kpc for distance of the 
Large and Small Magellanic Clouds respectively (Table|6ll. 

3.6. Results with different potential 

Calculations made with the updated potential are labelled with 
an asterisk in the tables and give the following results. The re- 
tention fraction fboimd exhibits a significant increase, especially 
for models lA* and IC*, while for the remaining ones the en- 



hancement is less conspicuous. In all cases the fraction of NSs 
retained by the disk is only slightly increased (Tabled. 

In all models with the updated potential, the scale height of 
the population is lower due to the increased masses of the disk 
and bulge, the decrease in z 1/2 being of the order of 10 - 20 per- 
cent (see the half density half thickness in Table|4|. 

The higher number of NSs retained by the disk implies also 
higher values of the surface and volume densities, together with 
the projected density for LOSs in the Galactic plane. For the 
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Fig. 11. Density profile on the Galactic plane (b - 0, upper 
panel) and along a meridian (1 = 0, lower panel) for Nsrar =10^ 
- Model IC. Cut-off distances are 30 (solid line), 10 (dotted) and 
3 (dashed) kpc respectively. 

same reason, even relatively fast NSs can be found in the disk, 
thus increasing the mean velocities on NSs (see Tables|4]and|5]). 

4. DISCUSSION 

Our calculations show that the distribution of birth velocities is 
the main factor driving the dynamics of NSs in the MW. We ob- 
tain substantially different values of fhound among different birth 
velocity models, with the shape of the distribution (position of 
the peak, bimodality, etc.) also playing a role in determining the 
final fate of bound NSs. 

The highest escape fraction, ~ 0.3, are obtained with H05, 
F06E an d F06L distributions. Th is value is lower than those 
found by lLyne and LorirneJ (1 19941) (hereafter LL) and A02 and 
similar to that inferred by H05. This is probably due to the 
fact that the velocity distributions proposed by LL and A02 
were obtained adopting the distribution of free electrons of 
iTavlor and CordesI Jl993h. Both H05 and F0 6 adopted instead 
the revised model of lCordes and Laziol (|2002|) for free electrons, 
which reduced estimates for distance, and thus velocity, of young 
pulsars. We performed a run with the distribution of A02, obtain- 
ing fhound ~ 0.54, confirming their results on the escape fraction. 

The adoption of the updated potential (higher mass of the 
Galaxy) implies higher escape velocities and hence only the 
fastest NSs, ~ 10-15 percent, can definitely escape from the 
MW. 

Albeit more than 70 percent of the NSs born in the MW are 
in bound orbits, the present-day number of NSs in the disk is 



/ 


^\ 




/■ 

°p-' - " 




/ 







■Ha^^^^^^UHBr- f 


\^ 


la": - ■ - •-- — ■ — "^^^^'^^^^^^^^^^ 




\.: ih..... ; 3 


90. .1.20. \ 




j 

la": - ■ - •-- — -— ^^^^^^^ 



Fig. 12. Sky maps of the projected density {Nsiar = 10*^) - Model 
IC. The cut-off distances are 30 kpc (upper panel), 10 kpc (cen- 
tral panel) and 3 kpc (lower panel) respectively. The density 
scale is normalized to the maximum density at 30 kpc. 



only a small fraction of the total, < 0.20. The remaining ones 
are found in the halo where they spend most of their life. This 
is a striking result but was not totally unexpected because, given 
their high spatial velocities, our synthetic NSs leave the disk in 
a short timescale, ~ 1 - 10 Myr Another remarkable finding 
is that the ratio of young to old NSs in the disk is very low: for 
each neutron star detected as a young active source there should 
be still more than 100 old NSs hiding in the disk. 

Our simulated NSs are born with significantly higher veloc- 
ities with respect to what is found in other works. In spite of 
this, our results for the half density half thickness show no sig- 
nificant differences with previous studies (except for the F06P 
distribution). Also, the local spatial density of NSs falls between 
those found by BR, BM, Z95 and that of P90, i.e. approximately 
between 1 and 5 x 10"'^(A^.5,flr/10'')pc"-'. This means that the 
nearest neutron star lies within ~ 10 pc from the Sun. 

The mean velocity is higher by at least a factor ~ 2 with 
respect to, for example, that found by BR, BM and Z95. Low 
velocity NSs (v < 50km s"') in the disk are a tiny fraction, 
/so ~ 0.001. This fraction grows by a factor ~ 10 in the LSR 
where /5^^* is ~ 0.05. Again results obtained with the F06P dis- 
tribution show a rather different behavior, with /50 <K 0.01 in 
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the GRF while in the LRF roughly half of NSs in the disk are 
in the low velocity tail (Tables |4] and |5]l- However, the effective 
weight of the low velocity tail of the distribution of birth veloci- 
ties is still matter for debate. 

Most of NSs, both bound and unbound, run away from the 
Galactic plane in a short timescale and form a halo which ex- 
tends well beyond the virial radius of the MW. The phase-space 
distribution of halo NSs clearly shows a separation between 
bound and unbound NSs. Unbound NSs become dominant at 
r ~ 500 kpc. 

The results presented in this paper will enable us to revisit a 
number of problems concerning isolated old NSs, like the accre- 
tion luminosity and its observability, the strategies for observing 
very close NSs, say within 100 pc and the optical depth of NSs 
in the perspective of using gravitational lensing to probe the pop- 
ulation. 
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Appendix A: Coefficients of the fits. 

We give here the best fit parameters for the surface density of 
the disk (Table [ATTT i, the height distribution (Table |A!2] i and the 
cumulative velocity distributions in the disk, both in the GRF 
(Table|A3]l and in the LSR (Table|A31i. 
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Table A.l. Surface density of the disk. 



Model 


lA 


IB 


IC 


ID 


lli 


lA* 


IB* 


IC* 


ID* 


m* 


Oo 


6.09 


6.48 


6.51 


6.60 


6.54 


6.24 


6.57 


6.57 


6.66 


6.61 


[xlO-'] 


-2.54 


-2.58 


-3.08 


-2.79 


-1.74 


-2.47 


-2.41 


-2.84 


-2.68 


-1.83 


Cl2 

[xlO-2] 


1.40 


1.57 


2.37 


1.83 


0.32 


1.05 


1.20 


1.96 


1.79 


0.58 


[xlO-'*] 


-5.69 


-8.06 


-12.70 


-9.20 


-2.28 


-2.99 


-5.25 


-10.22 


-10.07 


0.13 


[xlO-5] 


0.90 


1.67 


2.59 


1.82 


-0.34 


0.22 


1.00 


2.11 


2.25 


0.34 


Error 
[%] 


4 


3 


3 


3 


2 


5 


3 


2 


3 


2 



Table A.2. Distribution of heights. 



Model 


lA 


IB 


IC 


ID 


IE 


lA* 


IB* 


IC* 


ID* 


IE* 


bo 


142.3 


87.5 


123.1 


92.9 


112.7 


122.7 


73.7 


115.5 


91.5 


107.4 


h 


1.12 


1.20 


1.16 


1.22 


1.22 


1.14 


1.25 


1.17 


1.24 


1.24 




-136.3 


-83.7 


-120.0 


-91.0 


-111.9 


-117.1 


-70.6 


-112.7 


-89.9 


-106.8 


Error 


20 


41 


50 


63 


48 


22 


46 


63 


65 


56 



[%] 



Table A.3. Cumulative velocity distribution in the disk. 



Model lA IB IC ID IE lA* IB* IC* ID* IE* 



vo 

[kms '] 


214.1 


207.4 


201.1 


200.1 


207.2 


244.5 


241.2 


232.5 


232.5 


238.6 


n 


3.98 


4.56 


4.05 


4.63 


8.09 


4.03 


4.84 


4.16 


4.85 


8.21 


m 


3.98 


4.56 


4.05 


4.63 


8.09 


4.03 


4.84 


4.16 


4.85 


8.21 


Error 
[%] 


72 


91 


91 


95 


100 


63 


87 


93 


97 


100 



Table A.4. Cumulative velocity distribution in the disk (LSR). 



Model 


lA 


IB 


IC 


ID 


IE 


lA* 


IB* 


IC* 


ID* 


IE* 


[kms '] 


162.4 


130.5 


169.4 


133.5 


57.1 


178.7 


139.7 


182.3 


141.9 


64.4 


n' 


3.35 


3.33 


2.77 


2.57 


1.92 


3.34 


3.35 


2.70 


2.50 


1.94 


m' 


2.35 


3.33 


2.76 


2.57 


1.92 


3.34 


3.35 


2.70 


2.50 


1.94 


Error 
[%] 


73 


70 


44 


18 


26 


67 


62 


43 


12 


26 



